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We study the phase transition between the Coulomb hquid and the columnar crystal in the 3D 
classical dimer model, which was found to be continuous in the 0(3) universality class [l[. In 
addition to nearest neighbor interactions which favor parallel dimers, further neighbor interactions 
are allowed in such a manner that the cubic symmetry of the original system remains intact. We 
show that the transition in the presence of weak additional, symmetry preserving interactions is 
first-order. However the universality class of the transition remains continuous when the additional 
interactions are weakly repulsive. In this way, we verify the existence of a multicritical point near 
the unperturbed transition and we identify a critical line of unconventional transitions between the 
Coulomb liquid phase and the 6— fold columnar phase. 

PACS numbers: xxxx 
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Phase transitions are fundamental in nature and as a 
consequence they are relevant to every branch of clas- 
sical or quantum physics. They also serve as an im- 
portant playground for conceptual development of new 
ideas. A powerful approach to phase transitions and crit- 
ical phenomena is the Landau-Ginzburg- Wilson (LGW) 
theory ^ which is based on the concept of a local order 
parameter which characterizes a phase where some sym- 
metry is spontaneously broken. Then, the corresponding 
free energy of the system in the region around the crit- 
ical point is expanded in powers of the order parameter 
and several properties can be studied using renormal- 
ization group methods for continuous phase transitions. 
Recently, Senthil et al. Q proposed that phase transi- 
tions between phases with different non-trivial symme- 
tries can be, in principle, continuous without requiring 
fine-tuning. In particular, in the context of quantum an- 
tiferromagnets, it is thought that such a continuous phase 
transition from the Neel state to a valence bond solid 
can be reahzed contradicting the LGW prediction of 
a first order transition and its recent observation in re- 
cent simulations , an intermediate disordered phase or 
a fine-tuned direct continuous transition. 

In particular, phase transitions in systems composed of 
classical dimers at close packing could serve as an impor- 
tant example to illustrate the idea of a non-LGW thermal 
transition fit]. In such lattice systems, links are covered 
by close-packed hard-core dimers. The constraint of hav- 
ing precisely one dimer emanating from each site of the 
lattice allows for plausible connections to be made with 
associated gauge theories 0, capturing the essence of 
non-LGW transitions. In the case of the three dimen- 
sional (3D) cubic classical dimer model, the development 
of a featureless liquid with short-range correlations is not 
allowed even at infinite temperatures because then the 
dimer correlations fall off algebraically with the distance. 



being dipolar, and gapless collective modes (photons) are 
present [8|. This so-called Coulomb phase has also been 
found in other three dimensional model s [p and is rele- 
vant to some pyrochlore compounds [l^, ITll |. In the pres- 
ence of local attractive interactions, strong evidence was 
put forth P, that there is a continuous transition from 
the Coulomb liquid phase to a 6-fold degenerate columnar 
phase in the 0(3) universality class, which could support 
a non - LGW theoretical description. However, the ques- 
tion still remains open whether the transition is highly 
fine-tuned or weakly first-order. Several theoretical sce- 
narios have been suggested for the nature of this tran- 
sition [l^, [13] , proposing that the transition has an ef- 
fective SU{2) symmetry, and some evidence is linked to 
a mapping on a 2D quantum problem [13]. It is argued 
that this symmetry is not broken, since high-order terms 
are required to violate the symmetry. Here, we study 
the same problem in the presence of weak, same symme- 
try but apparently higher order perturbations, seemingly 
irrelevant in the RG sense. 

The system that we are interested in P, |^ is a cubic 
lattice of volume V = where L is the linear system 
size, covered by hard-core dimers and interactions which 
favor parallel alignment of dimers and respect the cubic 
symmetry. The energy is determined by the number of 
parallel dimers in x, y and z directions at different dis- 
tances but same column E ~ E= + E\ i 
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where a can be =, ]], //, i.e. pair of parallel dimers along 
the X, y or z directions and correspondingly = x, y, z. 
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FIG. 1: Phase diagrams of the described models in 
the presence of the interaction J2 (and J3 — 0) (inset: 
J3, with J2 = 0). In both cases, the lines shown separate 
the Coulomb liquid phase from the columnar phase. The point 
where J2(J3) is zero is consistent with the result reported 
m Ref 111. The dashed lines denote first-order transitions 
whereas the solid signify continuous ones. The lines through 
the points are guide to the eye. 

By n^ii) we denote the absence or existence (0 or 1 re- 
spectively) of dimer on the link (r^, fi+em) of the lattice, 
where e„i is the unit vector of the lattice in each of the 
three directions (m = x, y, z). Note that the dimer orien- 
tation is perpendicular to the direction along which the 
interactions are considered. 

The common property of the perturbations J2 , Js, is 
that they have the same sign with the original interac- 
tion, and they respect the cubic symmetry. We choose 
both J2 and J3 to be attractive, since if repulsive, there 
is additional frustration, similar to that of the ANNNI 
mode l [isl . . the study of which is the subject of future 
work 'l7!|. Given the order parameter of the transition, 

ma(r) = ^(-l)''°[na(r) - ^^(r - Cq)] , (3) 

the two sets of ground-states, classified according to their 
order parameter expectation value, favored by J2 and J3 
respectively, are essenially different. 

In this paper, we show that (I) Small perturbations 
that keep the original cubic symmetry of the Hamilto- 
nian intact and have the same sign (attractive), drive the 
transition from the Coulomb liquid to the columnar solid 
first - order, (II) The first-order transition gets stronger 
as the strength of the perturbation increases, (III) Our 
study is consistent with having a continuous transition 
for the system without the interactions J2, J3, or when 
they are repulsive. 

We perform calculations with lattices of linear size 
L=12, 16, 20, 24 and periodic boundary conditions, us- 
ing the Wang-Landau Monte Carlo algorithm Is] which 
is efficient in detecting first order transitions. The Wang- 
Landau scheme determines the density of states g(E) in- 
dependently of the temperature, therefore there is a di- 
rect access to the Gibbs probability histograms. Indepen- 




FIG. 2: The order parameter m (inset) and its sus- 
ceptibility for J2/J1 = 0.05. The visible "shoulder" feature 
in the specific heat does not appear here. The peak of the 
susceptibility, which signifies that the high coupling phase is 
columnar, scales as as one expects for a first-order transi- 
tion. 
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FIG. 3: Speciflc heat of the model with J2 = 0.05Ji. 

The peak scales asymptotically as ~ L^, expected for a first- 
order transition [1^. The "shoulder" feature, observed in 
Ref. [l[ is also present here, even though much weaker and 
goes away as the interaction strength gets bigger. In the inset, 
the energy histograms develop a 2-peak structure for J3 = 
O.OSJi, characteristic of a first-order transition. 



dent random walks are performed in different, restricted 
ranges of energy, by flipping plaquettes occupied by near- 
est neighboring parallel dimers. 

Each time a state of energy E is visited, we update the 
corresponding density of states by multiplying the exist- 
ing value by a modification factor / > 1, g{E) — > g{E)f. 
The density of states is modified until the accumulated 
histogram H{E) is relatively flat and it converges to the 
true value with an accuracy proportion to log(/). By re- 
fining continuously the modification factor (/ — > //\/2) 
and repeating the prescribed procedure, we achieve 10~^ 
accuracy for the density of states. 

The phase diagrams of the 3D dimer model in the pres- 
ence of J2 or J3 are shown in Fig. [TJ The transition tem- 
peratures are specified by using the characteristic fea- 
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tures that the energy cumulant Vl — 1 



has at 



a first order transition fio']. More specifically, Vl has a 
characteristic dip at the first-order transition point and 
takes the value 2/3 at all other couplings. The position of 
this dip converges exponentially fast with the system size 
and extrapolation gives an accurate estimate of {Ji/T)^. 
Our results are consistent, in the absence of any addi- 
tional perturbations, with the finding by Alet et al. of 
Jic = 0.597 [1]. As the interactions grow though, the 
transition becomes first order, with a decreasing correla- 
tion length. The transition moves to smaller couplings, 
as expected since the additional interactions are attrac- 
tive. For the smallest coupling available J2/J1 = 0.05, 
the specific heat, as it is shown in Fig. [31 shows the big 
"shoulder" feature that was observed by Alet et al. in the 
parent system, but is remarkably weaker, and gradually 
vanishes as the interaction gets stronger. As it is shown 
in Fig. m the characteristic 2-peak structure for the en- 
ergy probability histograms, gets widened as the interac- 
tion grows and at the same time the specific heat peak 
gets larger, for the same lattice size. Also, by tracking 
the distance between the peaks of the energy histograms 
as a function of the interaction strength J2/J1 or J3/J1, 
we found that there is a consistency of our result with 
having a continuous transition in the absence of any ad- 
ditional interactions, as it was found by Alet et al. 
The distance between the peaks gets closer to zero, if ex- 
trapolated to zero coupling J2 or J3, as shown in Fig. [5l 

We also studied the case where the interactions J2 are 
repulsive. We find, by performing simulations, that for 
small J2/J1 = —0.1, —0.2, —0.4 there is no signature of 
double - peak structure in the histogram structure, show- 
ing no signs of a discontinuous transition (cf. Supple- 
mentary material). Moreover, the scaling behavior of 
the susceptibility and specific heat peaks, signify that 
the transition has exponents that are close to the ones 
found in Ref. jH, {a/v = 0.75-1.0 and ^i/v = 1.75-2.0, 
giving v = 0.5 — 0.53 if we use hyperscaling) even though 
our simulated system sizes cannot permit us to conclude 
whether the transition for Jij J\ < is in the same uni- 
versality class. 

In Ref. ij] it was suggested that the critical theory 
of the classical dimer problem on the cubic lattice with 
nearest neighbor interactions is SU(2) invariant based 
on a mapping to a transition between a bosonic super- 
fiuid and a Mott insulator at fractional filling. By re- 
taining the same dimer ground states in the presence 
of further neighbor interactions, namely the six-fold de- 
generate ground states, the SU(2) invariance of the crit- 
ical theory should stay intact. In contrast, when the 
cubic symmetry is suppressed by invoking a symmetry- 
breaking potential, it has been shown that the transition 
becomes first-order [17]. 

In order to investigate how the presence of J2 and J3 
interactions affect the nature of the transition, we need to 
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FIG. 4: First-Order transition signature. The distance 
between the two peaks in the energy histograms gets bigger, 
as the interactions (J2) gets larger, as it is shown for L = 16. 
In the inset, the specific heat develops a larger peak (for the 
same lattice size L = 16) as the interaction gets larger. A 
very similar picture holds also for the case of attractve J3 
interactions. 
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FIG. 5: Weakness of the first-order transition for 
small perturbations. The distance between the two energy 
probability peaks when they have equal height, as a function 
of the interaction strength, for the two types of interactions 
we studied. Clearly, our study is consistent with a continu- 
ous transition for the parent system, since the slope becomes 
sharper near (0,0), making a non-zero value at J2,3/Ji = 
improbable. The dashed lines are guides to the eye. 



show how the coefficient of the 4-th order term in the GL 
functional of this transition changes. Even though the 
GL functional for such transitions cannot be derived in 
an exact way, we can use a crude approximation in order 
to investigate the rough effects of additional interactions. 
With this in mind, we write firstly the partition function 
of the system with just Ji interactions present, in terms 
of Grassmann variables [13, [13 '■ 

^ = j T^'n'Dri^ exp ^ NUj-q.-ql-qj-q] 

+ \^'^i^ijLmvlvjVjVkvivivl) , (4) 

ijkl 
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where Ji = 2{e''^ — 1), Mij represents the coordination 
array of the lattice and M-^l^ {p = 1, 2, 3) is the interac- 
tion matrix which is non-zero and equal to 1 when the 
four sites i, j, fc, I are the end points of interacting, with 
strength Ji, dimers. 

By using standard methods, we can integrate the 
Grassmann variables at the expense of introducing two 
Hubbard-Stratonovich fields and then, by defining the 
conjugate densities to the fields we can have the Gibbs 
potential in terms of the densities. Then we write the 
order parameter of the columnar order as ^^(r) = 
(—lY°m + Too and we identify one of the two Hubbard- 
Stratonovich fields by that, while the other Hubbard- 
Stratonovich field is represented by a spatially uniform 
function n. The Gibbs free energy density in terms of 
these fields reads: 

F = sJiTO^ + 6 JiTOq + 24JiTOon^ 



-Zn[36n2(Ji(m + 2mo) + 1)(Ji(-to + 2mo) + 1)] 



2JlTOo + 1 



(5) 



which we minimize with respect to toq and n, leaving m 
as the only order parameter. We need to emphasize here 
that by minimizing with respect to toq and n, we treat 
the dimer constraint (which is inherent in the grassmann 
structure of the theory) in a mean-field fashion. Such 
an approach is crude, similar in spirit to slave boson ap- 
proaches for t- J models Q , but it captures basic features 
of the model, such as tendencies under application of dif- 
ferent interactions. The Coulomb-liquid correlations are 
replaced by isotropic liquid correlations in this case. 
The free energy at the saddle - point can be expanded 



self-consistently in powers of to (17|. Such a mean-field 
expansion shows a continuous transition at — 0.51 
which we identify as a crude mean-field limit of the tran- 
sition in the absence of Jr for r — 2 or 3. In order to 
include the additional interactions, we perform a clus- 
ter expansion (where the extra terms are proportional to 
(±TO + moY where 2 > 3) to calculate the Gibbs poten- 
tial in a perturbation expansion in powers of Jr and ulti- 
mately, it takes the form [i3| r( Ji, Jr)/A^ = Co( Ji, Jr) + 

C2{Jl,JrW + Ci{Ji,Jr)m^ + C(i{Jl,Jr)m^ + ... . We 

find that the prefactor C4(Ji, Jr) is positive for values 
of Jr less than a critical value and becomes negative for 
larger values. The addition of more terms result in a 
renormalization of this critical value where the multicrit- 
ical point exists to lower values. The accuracy of the 
method is such that we cannot identify exactly Jr = 
as the multicritical point, nevertheless the change of the 
transition from second to first order can be explained. 

In conclusion, we investigated the phase transition be- 
tween a columnar crystal and a Coulomb liquid in a 3D 
dimer system. We showed that attractive perturbations 



which respect the symmetries of the original Hamiltonian 
drive the transition first - order. However, when the in- 
teractions are repulsive, the transition becomes continu- 
ous, seemingly in a different universality class than the 
one found in Ref. but a more extensive study is re- 
quired to characterize appropriately the nature of the 
transition. In this way, we might have given strong evi- 
dence for the nature of the non-compact CP^ universality 
class. 
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SUPPLEMENTARY MATERIAL FOR THE 
MANUSCRIPT: FIRST- ORDER VERSUS 
UNCONVENTIONAL PHASE TRANSITIONS IN 
THREE-DIMENSIONAL DIMER MODELS 

Here we present evidence that the transition in the 
repulsive case J2 < (keeping Ji > 0), is continuous, as 
it appears from the susceptibihty and specific heat data. 

The transition, up to the sizes that we have studied 
(maximum hnear size L — 32), appears to have similar 
characteristics with the transition found by Alet et ai, 
but the exponents seem to drift towards different val- 
ues. We cannot exclude the possibility that the tran- 
sition is actually in the same universality class, but up 
to the sizes we studied, such a conclusion is not obvi- 
ous. The transition has the same characteristics until J2 
becomes adequately negative (J2 ~ — 0.5Ji), where the 
6-fold columnar state becomes unstable towards the for- 
mation of 2 X 2 X 2 cubes filled with parallel dimers, mu- 
tually staggered. This phase is highly degenerate, since 
the cubes can fiuctuate between their different orienta- 
tions. So, for J2 < — O.5J1, a novel transition takes place 
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FIG. 6: Susceptibility of the order parameter x and 
the scaling of its maxima for J2/J1 — —0.4. The transi- 
tion is clearly continuous. 

between the Coulomb liquid at high temperatures and a 
highly degenerate state at low temperatures, character- 
ized by 2 X 2 X 2 cubes filled by parallel dimers, which 
can fiuctuate freely. However, this transition is going to 
be a subject of a future work. 

In Figs. [HI El we show the behavior of the system in 
the case J2/J1 = —0.4 which is well deep in the phase 
which extends for low J2, ranging from J2 — .1] to 
J2 ~ — 0.5Ji. The results on the critical exponents for 
smaller J2/J1 = —0.1,-0.2 are consistent with the re- 
sults we report in Figs. El [T] In Fig. |6] the specific heat 
shows an instability at Ji/T^ ~ 1.92 towards the 6-fold 
columnar liquid, and up to the sizes studied, the criti- 
cal exponent for the specific heat, C„ ^ cq + ciL"^" is 
compatible with a/iy ~ 0.75 — 1.0. Also, the suscepti- 
bility of the columnar order parameter is diverging with 
the system size according to x ^ L''/" ^ L^^'' where 
is the anomalous exponent of the order parameter. We 
find, as shown in Fig. [71 that j/v ~ 1.75 — 2.0. If we 
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FIG. 7: The specific heat and the scaling of its maxima 

for J2/J1 — —0.4. In this case a > 0, in contrast with the 
results of Ref. 01, but close to the values reported in Ref. 



use hyperscaling, namely that a = 2 — i/rf, wc find that 
v = 0.5 — 0.53. In conclusion, the transition in the repul- 
sive case is continuous and the exponents are close, even 
though seemingly not the same with the ones that are 
reported in [l| . The existence of a crossover towards the 
exponents reported in Ref. [2] might be a possibility, but 
the verification of this scenario requires a more extensive 
study in much bigger systems. 
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FIG. 8: Sample histogram for interaction J2 = O.OSJi 
and system size L — 24. For all energies below the maxi- 
mum possible energy, iSmax = Ji{l + J2/ Ji)L^ the histogram 
values are above 80% the average histogram value. 

The histograms, according to the Wang-Landau recipe, 
are designed to converge when they are 'fiat', closer than 
80% of the average histogram value for all recorded en- 
ergies (being relevant around the transition point) (cf. 
Fig. [8|). The procedure of making the histograms flat 
is repeated for 16 times, before finally trusting the form 
of the density of states (again, this is the typical proce- 
dure used in the applications of the Wang-Landau algo- 
rithm i^). 



[1] F. Alet, G. Misguich, V. Pasquier, R. Moessner and J.L. 

Jacobsen, Phys. Rev. Lett. 97, 030403 (2006). 
[2] D. Charrier, F. Alet and P. Pujol, Phys. Rev. Lett. 101, 

167205 (2008). 



[3] F. Wang and D. P. Landau Phys. Rev. Lett 86, 2050 
(2001). 



